Desynchronization in diluted neural networks 



O 

o 

(N 

in 



C/3 



Riidiger Zillmer,^'0 Roberto Livi,^'^'Q Antonio Politi,^^^'0 and Alessandro Torcini^- ^ •[! 

^ INFN Sez. Firenze, via Sansone, 1 - 1-50019 Sesto Fiorentino, Italy 
Dipartimento di Fisica, Universitd di Firenze, via Sansone, 1 - 1-50019 Sesto Fiorentino, Italy 
'^Seztone INFN, Unita' INFM e Centro Interdipartimentale per lo Studio delle Dinamiche Complesse, 
via Sansone, 1 - 1-50019 Sesto Fiorentino, Italy 
'*Istituto dei Sistemi Complessi, CNR, CNR, via Madonna del Piano 10, 1-50019 Sesto Fiorentino, Italy 
^Centro Interdipartimentale per lo Studio delle Dinamiche Complesse, 
via Sansone, 1 - 1-50019 Sesto Fiorentino, Italy 

The dynamical behaviour of a weakly diluted fully-inhibitory network of pulse-coupled spiking neu- 
rons is investigated. Upon increasing the coupling strength, a transition from regular to stochastic- 
like regime is observed. In the weak-coupling phase, a periodic dynamics is rapidly approached, 
with all neurons firing with the same rate and mutually phase-locked. The strong-coupling phase is 
characterized by an irregular pattern, even though the maximum Lyapunov exponent is negative. 
The paradox is solved by drawing an analogy with the phenomenon of "stable chaos", i.e. by ob- 
serving that the stochastic-like behaviour is "limited" to a an exponentially long (with the system 
size) transient. Remarkably, the transient dynamics turns out to be stationary. 

PACS numbers: 87.19.La,84.35.+i,05.45.Xt 
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I. INTRODUCTION 



During the last few years it has become increasingly 
clear that understanding the behaviour of many different 
systems passes through the comprehension of the dynam- 
ics of complex networks Ij. This is, for instance, the case 
of metaboUc systems, genetic networks, the immune re- 
sponse system and neurobiological structures . A par- 
ticular challenge is represented by the need of unravelling 
the mutual connections between network structure and 
dynamical properties. Very little is in fact known about 
the expected classes of behaviour and their stability prop- 
erties even in systems of globally coupled identical oscilla- 
tors. For this reason and under the assumption of a struc- 
tural stability of the possible scenarios, it is, therefore, 
instructive to investigate simple models, such as diluted 
neural networks of pulse-coupled neurons. Since it seems 
that inhibition plays a major role in determining the dy- 
namics of single neo-cortical piramidal neuron ^ as well 
as of cortical networks we have chosen to examine 
a network of inhibitory coupled leaky integrate-and-fire 
neurons. More precisely, we consider the model proposed 
in Ref. 4], where Jin, under fairly general conditions, in- 
vestigated analytically the convergence towards a peri- 
odic pattern. In this paper we study the diluted version 
of this model, showing that even though the dynamics 
is characterized by a negative maximum Lyapunov ex- 
ponent 0, irregular and exponentially-long transients 
are typically observed for a sufficiently strong coupling 
strength. This phenomenon is somehow analogous to 
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what observed in diluted networks of pulse-coupled os- 
cillators with delay , although the transient therein re- 
ported are chaotic in the typical sense of the word. More 
precisely, we find that for small coupling amplitudes the 
dynamics converges, after a short transient, towards a 
synchronized state with all neurons firing with the same 
rate, but with their phases approximately uniformly dis- 
tributed. This can be understood by first referring to 
a homogeneous network of globally coupled neurons. In 
that context, the mean field is found to induce an effective 
repulsive interaction between the neurons; as a result the 
asymptotic regime is characterized by an evenly spaced 
sequence of spikes which is known in the literature as a 
splay state l7|. As a result of random dilution, one can 
imagine that inhomogeneities in the mutual interactions 
arise which in turn lead to small nonuniformitics in the 
interspike intervals. 

On the other hand, for sufficiently large coupling am- 
plitudes, stochastic- like transients are observed, whose 
duration is exponentially long with the network size. 
Since various indicators show that this regime is station- 
ary, it is logical to conclude that in infinitely large net- 
works it represents a perfectly legitimate thermodynamic 
phase. This is analogous to the active phase in directed 
percolation, whose life-time is finite in finite systems. 
Even more stringent is the analogy with "stable chaos" , 
a kind of irregular behaviour discovered in coupled map 
lattices!^ and characterized by negative Lyapunov ex- 
ponents. Since it is believed that such a pseudo-chaotic 
behaviour is sustained by discontinuities of the mapping 
ruleQ, it is natural to expect this to be true also in the 
present case. Quite consistently, we observe that in the 
presence of disorder, where neurons are no longer equiv- 
alent to one another, changes in the firing order are ac- 
companied by discontinuities in the evolution rule. 

Altogether, the transition manifests itself as a col- 
lective de-synchronization phenomenon. Unfortunately, 
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a "microscopic" linear-stability analysis does not allow 
identifying the threshold, since all trajectories are asymp- 
totically stable in both regimes. Moreover, direct numer- 
ical simulations are not very effective either, due to the 
difficulty of simulating large networks, so that the study 
of the critical behaviour is an even more difficult task. 
However, the evidence of two distinct phases is rather 
convincing and, furthermore, the introduction of a suit- 
able space-time representation suggests a true analogy 
with directed percolation that will be worth exploring in 
more detail. 

In Sec. m we introduce the model and the variables, 
while the homogeneous case of fully coupled networks is 
analytically investigated in Sec. IIIII where we also intro- 
duce a one-dimensional description of the dynamics that 
applies exactly in the thermodynamic limit. In Sec. IIVI 
the transient dynamics of weakly diluted networks is dis- 
cussed, both by determining the transient length and 
analysing its stationarity properties. The resulting two- 
phase scenario is then summarized in Sec. El where we 
also present an effective stochastic description. Further 
comments about open questions and indications for fu- 
ture studies are finally presented in the conclusions. 



II. THE MODEL 

In this paper we investigate a system of N leaky 
integrate-and-fire (LIF) neurons, analogous to the model 
discussed in Ref. Q . The state of the i-th neuron is fully 
determined by the membrane potential Vi{i) and obeys 
the differential equation 



5™^), (1) 



where t is the membrane time constant, C is the 
suprathreshold input current (referred to a unitary mem- 
brane resistance), and W is the reversal potential. When- 
ever the potential Vj{t) reaches the threshold value 8, it 
is reset to i? < 8, and a spike is sent to and instanta- 
neously received by all connected neurons at time tj™'' 
(the superscript m enumerates the firing events of the 
j-th neuron). The net result of a received spike is that 
the membrane potential of the i-th receiving neuron is 
decreased according to the transformation 



Vl + W ^{V, + W)eM-9^J)■ 
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As a consequence of the inhibitory connections, the po- 
tential VI can go below the reset value R, but Eq. Q 
shows that — Wisa, true lower bound. However, for small 
coupling values and in the absence of clustering phenom- 
ena (as in the present manuscript), it turns out that Vi 
is essentially bounded between R and 8. The last ingre- 
dient defining the system dynamics is the connectivity 
matrix gij. Following the recent literature on randomly 



connected directed networks the coupling strength 

is scaled to the connectivity of the receiving neuron. 



9ij 



G I li , if i and j are coupled 
, otherwise 



(3) 



where G is the coupling constant and li is the number 
of incoming links to neuron i . In other words, we con- 
sider the simplest type of disorder, determined just by 
the presence/absence of links between neurons (notice 
that self-interactions are excluded, i.e. ga = 0). As we 
are interested in studying networks with a given frac- 
tion Vra of missing links, there are in principle different 
ways of doing that: (i) each link is cut with a proba- 
bility r™; {ii) the total number of cut links is fixed 
deterministically [Nm — rmN{N — 1)); (iii) the number 
of cut links per neuron is fixed deterministically. Since 
preliminary simulations performed according to the three 
philosophies have given qualitatively similar results, we 
have eventually decided to restrict our quantitative anal- 
yses to the second approach. Moreover, since we aim at 
understanding the possibly qualitative changes induced 
by disorder in the neural network dynamics, we limit our- 
selves to analysing the case of weak disorder. This is done 
by choosing small values for the fraction r„j (typically 

Tm = 0.05). 

Most of the parameter values are set according to the 
current literature (see e.g., 0,^3), namely: t = 20 msec, 
G = -45mV, W = 63mV, 8 = -52mV, and R = 
—59 mV. As a matter of fact, the only free parameters of 
the model arc the fraction and the coupling constant 
G, which tunes the strength of the inhibitory coupling. 

The dynamical equation ^ can be recasted into a sim- 
pler form by introducing the dimensionless parameters 

t^i/T,v, = (F, - i?)/(e - i?) , 

c = (C - i?)/(8 - R) , w ^ (W + R)/{e - R) . 

This amounts to rescaling 8 and R to the values 1 and 0, 
respectively. Moreover, for the above choice of parameter 
values, c = 2 and w = 4/7. 
As a result, Eq. ^ reads 



c-Vi-{vi + w) ^ ^ gij S{t- ij"-') 



(4) 



Since w > 0, the coupling is fully inhibitory, i.e., an in- 
coming spike lowers the potential of the receiving neuron 
and hence the coupling inhibits firing. 

At variance with the original model , where both the 
threshold currents Gi and the coupling constants gij are 
assumed to be randomly distributed, here the only source 
of disorder is the presence/absence of inhibitory connec- 
tions. This choice is dictated by our interest in relatively 
simple structures to better understand the possible sce- 
narios. However, the most important difference concerns 
the coupling constant: no dependence on the system size 
is postulated in Ref. 4] , while an inverse proportionality 
to the number of incoming connections is assumed here. 
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This latter normalization is more suited to investigate 
the large iV-limit, since it guarantees that the coupling 
strength remains comparable to the amplitude of the 
force field {c — v). 

By following the approach of Rcf. 4] , we map the orig- 
inal model onto a discrete-time map. Let Vi{n) denote 
the membrane potential of the i-th neuron immediately 
after the n-th spike. Until the next spike is emitted, the 
network evolution corresponds to an exponential relax- 
ation of each Vi towards c. One can easily infer from Q 
that the time interval Ati needed by the z-th neuron to 
reach the firing threshold = 1 is 

At,(7i) :=lnr,(n) with r,(n) = ^— . (5) 

c — 1 

Let now k{n) denote the label of the neuron characterized 
by the shortest time, i.e.. 



^k{n) = niin{ri} . 



(6) 



By now including the effect of the next spike, the network 
dynamics can be be transformed into a discrete map for 
the quantities ri(n). In the large N (and, accordingly, 
£i) limit, one obtains 



c-1 



(7a) 



ri{n + 1) =: gika+ {1 - gik)——— for i 7^ fc , (7b) 

1 k[n) 



where 



c + w 

c-1 



(8) 



The process can be iterated by finding the minimum 
among the Ti(n + 1) and so on. Therefore, Aij,(„-)(n) 
represents the time interval between the n-th and the 
(n -I- l)-st spike; following the standard notation, it will 
be denoted with iisi(«)- On the other hand T(m) de- 
notes the time elapsed between the (m -I- l)-st spike and 
the previous spike emitted by the same neuron. 



III. THE HOMOGENEOUS CASE 

In this section we discuss the dynamics of homogeneous 
networks, i.e., gij = G/£i (for i ^ j) and £i — N. In 
the absence of any disorder (globally coupled identical 
neurons), any initial ordering of the neuron spikes will 
persist at all times, because it cannot be changed by the 
dynamical rule. 

Accordingly, after a short transient, the trajectory con- 
verges towards a periodic orbit characterized by equis- 
paced spikes, tisi(n) = T/N and T{m) = T for all n, m. 
By formally interpreting the separation between the spike 

time t^™^ of the jth neuron and the preceding spike-time 
) of some reference neuron as a phase variable, one can 



summarize the scenario by stating that the single-neuron 
phases repel each other until an equilibrium state, charac- 
terized by a maximal and uniform separation, is attained. 
In the literature, this alignment is called splay state and 
it has been found in various models of globally coupled 
oscillators . Such a state can be viewed as the opposite 
situation of a fully synchronized state, where all oscilla- 
tors share the same phase. It is curious that this regime 
is attained in the presence of inhibitory coupling, since 
this property is usually associated with the propensity 
to synchronize instead. In fact, splay states have been 
found in globally coupled LIF neurons in the presence of 
excitatory coupling with spikes of finite width [T^ . 



A. Stationary solution 

Since the ordering of the neuron potentials, Vi, does not 
change in time, one is free to label the spiking neuron in 
such a way that k{n -I- 1) = k{n) + 1 mod N. It is also 
convenient to choose a "moving" frame, namely j = i — n, 
since the label of the spiking neuron remains constant and 
can thereby be set equal to 1, without loss of generality. 
The resulting map reads 



a G , 

r,_i(n + l) = — +(1- 



N J Ti{n) ' 



(9) 



with the boundary condition Fjv = c := c/(c— 1). This 
equation admits a stationary solution, which corresponds 
to a regime of evenly spaced spikes, i.e. a splay state 0. 

The fixed point of © is obtained by solving the recur- 
sive relation 



f _aG 



1 - 



T + G 
N 



r 
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where Fi is approximated with 1 + T/N . In fact, from 
the very definition, one has InFi = ijsi ~ T/N, so that 
T is the (constant) single neuron interspike interval. By 
iterating backward the above equation from the initial 
condition Fat = c, one obtains 

aG 



N-n — 



T + G 



1 _ ^-{T+G)n/N 



AT+G)n/N 



(10) 



The value of T can be finally determined by imposing the 
condition 

fi^l + T/N . (11) 
By neglecting first order corrections in 1/N , one obtains, 



1 



aG 
T + G 



1 



-{T+G) 



+ ce 



-{T+G) 



(12) 



B. Linear Stability 



Let us consider the evolution of an infinitesimal pertur- 
bation 7j (n) of the fixed point T j . From the lineariza- 
tion of Eq. @ it follows that 7^ [n) satisfies the recursive 
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relation 



7j_i(n+l) 



T + G 



7jH 



r 

TV 



rj7i(n) 



(13) 

where higher order corrections in are neglected. Be- 
cause of the reset F^v — c/{c — 1), it is natural to im- 
pose the boundary condition — . By further setting 
jj{n + 1) = fi"fj{n), one obtains the eigenvalue equation, 



(l-§)r,7i-/i7,-i=0, (14) 



where 



1 - {T + G)/N 



By iterating the above equation and imposing the bound- 
ary condition, we find that jl satisfies the equation 



T 



N-2 



3=0 



By inserting the expression of i^N-j (see Eq. (|1U|| ') and 
using the condition H12|) . we find that for 1 



Inc 



(15) 



which yields the expression of the Lyapunov exponent of 
the discrete map 



A = ln|M| = -^[r + G-ln5] 



(16) 



The Lyapunov exponent of the original system is finally 
obtained by rescaling time with the intcrspike interval 
T/N, 



A= ^ln|Ai| = -1 jr- 



(17) 



In Fig. n we plot A and compare it with numerical 
results for different values of the coupling G . It is worth 
recalling that T = T{G) as from Eq. (P|l . 



C. Continuous approach 

In the large N limit, the discrete time (n) and "space" 
(j) variables can be transformed into continuous ones by 
introducing 

j/N , X e [0,1] , r = njN . 
As a result, by neglecting 1/iV^ terms, Eq. (jSJ writes as 

dv or 



dx dr 



-{T + G)T + aG ■ 



(18) 




FIG. 1: The maximal Lyapunov exponent A corresponding 
to for Tp = 0, c = 2, w = 4/7 versus the coupling G. The 
dashed line refers to the analytical expression 11711 and the 
symbols to numerical estimates. 



where we have made use of Eq. (|11|) . The boundary con- 
ditions read now r(0, r) = 1 and r(l,r) = c. The sta- 
tionary solution corresponds to the fixed point of and 
reads 



f(x) 



1 



aG 
T + G 



AT+G)x 



aG 
T + G 



(19) 



The period T is defined by imposing the boundary con- 
dition at a; = 1, i.e., 



aG 

T + G 



'.G 



T + G 



This result is consistent with Eq. (|12|l . Thus, we see 
that the evolution equation of the globally coupled ho- 
mogeneous network can be reduced to a 1-1-1 dimensional 
partial differential equation. We shall see later that this 
analogy proves fruitful to interpret the phase transition 
discussed in the following sections. 



IV. TRANSIENT DYNAMICS 

From now on we study the dynamics of model Q when 
a small fraction of links is removed. It turns out that 
in this case a much more interesting and complex dynam- 
ics may emerge when the coupling strength is increased, 
even though the Lyapunov exponent remains negative. 
However, at variance with the previous section, here we 
can rely on numerical simulations only. 

We expect that the weak amount of disorder induced 
by the random pruning of directed links reduces the cou- 
pling strength among the neurons and correspondingly 
increases the value of A, with respect to the fully cou- 
pled case. The results plotted in Fig. ^ show that this 
is indeed the the Lyapunov exponent is found 

to increase with the pruning ratio . Nevertheless, A 
remains negative up to = 0.2. Moreover, numerical 
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simulations indicate that A remains finite in the thermo- 
dynamic hmit N ^ oo. 



o-oA'=100 
<>■-<> A'=200 
A.Ayv=300 




0.15 0.2 



FIG. 2: The maximal Lyapunov exponent A for different val- 
ues of the cut ratio rm and for three values of the system size. 
The data refer to a coupling G = 2. 

As a result, sooner or later the dynamics must con- 
verge towards a stable periodic orbit. In Ref. 4], Jin has 
derived an upper bound for the transient length in the 
case of generic disorder. More precisely, he finds that the 
number of spikes preceding a periodic pattern is 



P « 7V« + - 1) 



(20) 



where, in the large N limit, and with reference to our 
notations, 



q = -\n (A/8a^) /g,nin + 1 , 

gmin = min{.gij} , and finally, 

I'd 

A= min min [rj(n) - rfe(„)(n) 

n=l,...,oo j^k(n) 



Although the N dependence expressed by Eq. super- 
ficially looks like a power law, it is in fact much faster. 
First of all, if the coupling strength scales as (as it 
is assumed here), the exponent q grows linearly with N. 
Incidentally, this is true also in the fully coupled regime, 
where we have seen that the transient is, in reality, much 
shorter. A second reason is that in systems with some 
pruning (like here), gmin is equal to and the exponent 
q would be actually equal to infinity. Finally, A, which is 
the minimal distance between the membrane potential of 
the two neurons that are closest to the firing threshold, 
cannot be larger than This leads to an additional 

logarithmic growth of the exponent q with N. Accord- 
ingly, we conclude that the analytic estimate contained 
in Ref. applied to our setup is much too large an over- 
estimation of the transient length to help understanding 
when and whether a phase transition can occur. 

The transient duration tn is determined as the smallest 
time for which the actual configuration of the membrane 
potential is e-close to a previous state. 



where 1 < i < TV and 1 < T < t - 1. The choice of the 
threshold is not crucial, since after the maximal distance 
is on the order of 1/N, it starts converging exponentially 
fast to 0. Accordingly, by this procedure, we do not only 
determine the transient length but also the periodicity T 
of the asymptotic solution. 

The dynamics (0J is usually simulated by starting from 
a random initial condition with the values of the mem- 
brane potentials Vi{0) uniformly distributed in the inter- 
val [—w, 1]. In the following we present a detailed analysis 
carried out for networks with = 0.05. 

In order to deal with a more reliable quantity, we com- 
pute the average length of the transient (ttr) (here and 
in the following, angular brackets denote an average over 
both realizations of the network and different initial con- 
ditions of the membrane potentials). In Fig. 13 we have 
plotted the average transient versus G for two different 
values of the system size {N = 100 and 200). For G < 1, 
(iti) slightly decreases for increasing coupling strength. 
This is simular to what observed in the fully coupled net- 
work over the entire range of variation of G "4^. On the 
other hand, by further increasing the coupling strength 
(tti) exhibits a sudden increase. More important is that 
the rate of increase is significantly larger when TV doubles. 
Altogether this data suggests the existence of two dis- 
tinct phases approximately corresponding to G smaller 
and larger than 1. 
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niin{t| max \vi{t + T) — Vi{t)\ < £ , } 



FIG. 3: Average transient length (ttr) vs. G for two system 
sizes. 

A more precise characterization of the dynamics is ob- 
tained by investigating the scaling behavior with iV. In 
Fig. 0^, we see that for G = 0.5, the transient length 
increases linearly with the system size TV, while the av- 
erage period (T) of the asymptotic attractor remains al- 
most constant. At the same time, in Fig. QJ), we see 
that for G = 1.8 and G = 2.5, the transient grows ex- 
ponentially fast. Finally, it turns out that for G = 2.5 
also the average period grows exponentially. We consider 
this as a preliminary indication that at larger coupling 
strength another transition may occur. However, here 
we focus our attention on the qualitative changes occur- 
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ring for G ~ 1 where only the transient starts growing 
exponentially. 




FIG. 4: Average length {tti) (solid line) of transient and corre- 
sponding period (T) of the periodic attractor (dashed line) vs. 
system size: (a) G = 0.5 ; (b) G = 1.8 (circles) and G = 2.5 
(diamonds). Note in (b) the logarithmic scale of the vertical 



Another important point to be carefully investigated 
are the statistical properties of the transient dynamics. 
In particular, we start from the stationarity that can be 
analysed by computing the the so-called coefhcient of 
variation (CV), that is determined by subdividing a time 
series x{n) , n=l,2,...,L, into a sequence of windows 
of constant duration I, and then computing 

ml 

x{m)^- ^ a;(i) , m = 1,2, . . . , 



a{ra) = 



i=(m-l)i+l 
ml 

I 



^ ml 

7 (^(*) ^ x{m)f 



i=(m-l)i + l 



1/2 



The coefficient of variation is thereby defined as 
CV{m) ~ a{m)/x{m) . 



(21) 



The parameter I must be chosen in such a way that each 
window contains a sufficient number of samples to reli- 
ably compute the local average and the standard devia- 
tion. 

In our case, the time series of interest is the sequence of 
single neuron interspike intervals T, obtained by iterating 
map Q. All intervals are recorded over a time window 
of length I — ION, so that each neuron fires on average 
ten times. In order to reduce statistical fluctuations, the 
CV has been further averaged over 100 different initial 
conditions for a fixed network realization. In Fig. |31 we 
have plotted the time evolution of {CV{n))^ for different 



values of G. Moreover, for each value of G, the behav- 
ior of (Cy(n))j is presented for two different network 
realizations to convince the reader that a further aver- 
aging on the disorder is not truly required to conclude 
about the stationarity of the regime. Indeed, the main 
point is that for G < 1 the transient is nonstationary, as 
{CV{n))^ vanishes exponentially while approaching the 
periodic attractor. On the other hand, for G > 1.2, 
(Cy(ri))j approaches a finite value, thus confirming that 
the "transient" regime is stationary and one can thereby 
meaningfully speak of an invariant measure and pose the 
question of determining its correlation properties. 
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FIG. 5: Solid lines: Coefficient of variation as function of 
integer time for A'^ — 1000 and (in ascending order starting 
from the bottom) G = 0.9, 1.0, 1.1, 1.2, 1.3 . Stars: The same 
for a different disorder realization and G — 1.0, 1.1, 1.2 . The 
straight line corresponds to the saturation value at G — 1.2 . 

We have also performed other stationarity tests based, 
e.g., on the analysis of distributions of temporal distances 
of neighboring points in state space (e.g., see Q). They 
confirm the stationarity of the transients above the tran- 
sition and confirm the existence of a transition point in 
the interval 1.2 < G < 1.3. As these tools do not pro- 
vide further insight we do not further comment on their 
outcome. 

The transient length strongly depends on the initial 
condition. In Fig. we show the probability distribution 
function P(itr) for a given realization of the disorder in 
both phases. For small G, P(itr) can be confidently fitted 
with an inverse Gaussian (see the inset of Fig.EJ, which is 
the typical statistical distribution of biased escape-time 
problems 'iSl]. This confirms once more the nonstation- 
arity of the corresponding transient dynamics, since it 
is driven by a systematic drift towards the asymptotic 
periodic state. Above the transition, P(ttr) is Poisso- 
nian: the transient statistics can be therefore interpreted 
as a typical escape process through an activation barrier 
from a locally equilibrated system. Again, this evidence 
supports the conjectured stationarity of the transient dy- 
namics for G > 1 (in the large TV limit) . 

Finally, P(itr) does not change qualitatively for differ- 
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FIG. 6: Normalized histogram (random initial conditions at 
fixed disorder, A'^ = 65) of transient lengths for G = 2.1 . The 
solid curve is an exponential, the dashed curve an inverse 
Gaussian with equal mean and variance. The inset shows the 
distribution with an inverse Gaussian fit for G = 1 . 

ent realizations of disorder, although sizeable quantita- 
tive differences can be found for small values of N. How- 
ever, our simulations indicate that the probability dis- 
tribution P{ttr) is a self-averaging quantity since, upon 
increasing iV, sample-to-sample fluctuations become in- 
creasingly small. This is absolutely clear for G ~< 1, 
when it is possible to compute the transient time even for 
moderately large values of 0(10"*). Above G — 1, 
where direct simulations are almost undoable, there is, 
however, no reason to expect a qualitative change, given 
all the evidence of an ergodic dynamics. 

Altogether, the "transient" dynamics observed for 
G «> 1: (i) is characterized by a negative Lyapunov 
exponent; (ii) is effectively stationary; (iii) lasts for ex- 
ponentially long (with the system size) times. These are 
the distinguishing properties of "stable chaos", a phe- 
nomenon extensively investigated in coupled map lattices 
[i| , where similar features of the transient statistics have 
been observed However, it is worth recalling that 

stable chaos has been observed also in chains of forced 
oscillators [l^ and recently uncovered in the ID hard 
point gas of diatomic particles . 

In all such instances, stable chaos is associated with 
the presence of discontinuities in phase space. In the 
coupled map models, the discontinuities are transparent 
in the definition of the local piece- wise linear maps 0. 
In the hard-point gas they arise in connection to three- 
body collisions. Around any configuration leading to 
one such multiple collision, the ordering of the two-body 
collisions changes abruptly. The non-commutativity of 
the collision itself induces therefore a discontinuity in 
the dynamics that is conceptually similar to those pos- 
tulated in the coupled-map lattices. The same hap- 
pens in the present context. In fact, let us con- 
sider the time-evolved one-parameter family of configura- 
tions, v^{t) = {<(i), vl^{t), <(0, <+i • • ■ , V^^it)} 



where Vj (0) = uj + vSjk and 5jk is the Kroenecker delta 
function. It is easy to convince oneself that w^^ (t) breaks 
into two disconnected parts when v'j^{t) = 1 if is such 
that f^(i) = ffe_|_i(i) and there is only one connection 
between neurons k and fc -I- 1. The discontinuity is due 
to the fact that only one of them inhibits the other. It 
is therefore interesting to verify that, in agreement with 
the past observations, the presence of discontinuities in 
the phase-space is a condition for the onset of a "stable 
chaos" dynamics. If no anomalous transients arise for 
small coupling strengths, it is because the condition is 
necessary by not at all sufficient. 



V. TWO DYNAMICAL PHASES 

In this section we study the properties of the dynam- 
ical phases and perform a preliminary analysis of the 
transition. For G < 1, the evolution is similar to that 
of the fully coupled case. After a short transient time 
the system converges to a state characterized by a se- 
quence of N spikes (emitted by the N neurons), which 
periodically repeats itself (see Fig. [7J|. All neurons fire 
with the same pace {T{m) = const.), but their phases 
are no longer equispaced (i.e. tisi{n) varies in time). In 
other words, the asymptotic solution is a phase-locked, 
i.e. synchronized, state. In the following, this dynamical 
regime will be denoted Locked Phase (LP). The main dif- 
ference with the fully coupled case is that different spike 
sequences are not equivalent to one another. In the limit 
of identically-coupled neurons the dynamics is invariant 
under any permutation, but this is no longer true as soon 
as some degree of heterogeneity is introduced in the net- 
work connectivity. As a result, for rm ^ there exists an 
exponentially large number of different periodic attrac- 
tors. More precisely, the total number of attractors can 

be estimated to be of the order of A^!/(nd=i -^dO; where 
Nd is the number of neurons connected to the same (in- 
coming and outgoing) neurons in the network and M is 
the number of these equivalence classes present in a given 
network realization. 

For larger coupling strengths, an irregular dynamical 
regime arises that we call Unlocked Phase (UP), because 
the mutual ordering keeps changing during the entire 
transient, as one can, e.g., see in Fig. |S1 where slow but 
systematic pattern adjustments can be seen in the whole 
time range. 

One can better elucidate this dynamical phase by look- 
ing at the relative potential differences. Without prej- 
udice of generality, we choose one of the N neurons, 
say neuron i, and consider the subsequence of its fir- 
ing events. Given the spike emitted at time t, we com- 
pute the difference between the last ISI, r*^*\ of that 
neuron and the average ISI of the j neurons firing im- 
mediately before and after time t. We denote this time- 
dependent quantity by A0(*^(f). The value of j has to 
be chosen neither too large, to avoid inclusion of spikes 
from the same neurons, nor too small, to get rid of lo- 
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FIG. 7: Firing pattern (index of the firing neuron vs. time) 
of a typical periodic attractor for G = 2, A*' = 50 in the LP; 
k is the index of the firing neuron. 
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FIG. 8: Firing pattern for G = 2, iV = 50 in the UP; k is the 
index of the firing neuron. 



cal fluctuations. We have checked that the choice j — 5 
is a reasonable compromise. The integrated phase-shift 
A$(') {t) = /g A0(^) {t') is plotted in Fig.Elfor a generic 
initial condition and a typical realization of the disorder. 
In the LP, after a short transient, A<i>'^*-' collapses onto 
a constant value, which indicates that the system has 
converged to a periodic firing sequence with the neurons 
firing with the same pace. On the other hand, in the UP, 
A$('^ wanders erratically, performing a sort of unbiased 
random-walk. We can indeed confirm that at least for 
G < 1.8 and within statistical fluctuations all neurons 
are characterized by the same average firing rate. 

A spontaneous symmetry breaking induced by disorder 
appears at larger coupling strengths when we can, e.g., 
find periodic orbits with different neurons exhibiting dif- 
ferent firing rates (see Fig. 1101) . This phenomenon is, 
however, not connected with the onset of exponentially 
long transients that appears for significantly smaller G 
values. 

Further information about the two dynamical phases 
can be gained from the introduction of a suitable space- 
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FIG. 9: Time evolution of the local phase difference A<E>'''(f) 
for G = 1.0 (solid line) and G = 1.5 (dashed line). We 
consider a system with A'' = 1000 and diluition = 0.05 . 
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FIG. 10: Firing pattern associated with a periodic attractor 
in the UP for G = 2.0, iV = 50 . The period is indicated by 
the vertical lines. 



time representation. Having verified that all neurons sta- 
tistically behave in the same manner, we can split the 
entire spike series into sequences of N consecutive events 
and label the sequences with an index m, which naturally 
plays the role of time. On the other hand, the index i, 
labelling the position within a single sequence, plays the 
role of space. Finally, we introduce a binary variable 
6(to, i) to distinguish two cases: whenever the neuron 
emitting the z-th spike in the m-th sequence is the same 
neuron which emits the i-th spike in the (rn — l)-st se- 
quence, we set b{m,i) — 0; otherwise h{m,i) = 1. The 
resulting patterns are reported in Fig.^] where 0, 1 are 
coded as black and white colours, respectively. In the 
UP, the pattern is very irregular and no order seems to 
emerge even after a very long time lapse. Conversely, 
in the LP, black tends to prevail quite soon, indicating 
that the system rapidly converges to a fixed firing se- 
quence, i.e., to a periodic attractor. The most interest- 
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ing features of this representation are the white defects 
propagating in the black background. They indicate that 
the firing rate of some neurons is temporarily slower or 
faster than that of the "neighbouring" ones. These pat- 
terns are strongly reminiscent of a Directed Percolation 
(DP) transition. Actually, defects tend to be eliminated 
in the LP, which is analogous to the absorbing phase of 
DP. On its turn, the UP shows similarities with the so 
called "active" phase of DP, ruled by a persistent defect 
dynamics. It is remarkable that the dynamics of an al- 
most globally coupled system exhibits distinctive features 
of one-dimensional systems with short-range interactions, 
as it is the case of 1-1-1 contact processes. The analogy ap- 
pears less astonishing after recalling that the evolution of 
the fully coupled model is described by a suitable partial 
differential equation (see Eq. (|18|l ). However, including 
the role of disorder and proving a possible analogy with 
directed percolation is not an easy task. We prefer to 
leave the elucidation of this problem to future analyses 
and here we concentrate our efforts to obtaining a refined 
characterization of the two dynamical phases. 

In particular, we have determined the distribution and 
the time autocorrelation function of the single-neuron in- 
terspike time interval T^*) (m) . In the UP, the normalized 
autocorrelation function Cisi is plotted in Fig. ^1 It ex- 
hibits an exponentially fast decay, typical of both chaotic 
and stochastic regimes: memory is lost after typically 3 
to 4 firing events. This quantitative confirmation of the 
"irregularity" is yet another element strengthening the 
analogy with stable chaos. The "stochastic" -like charac- 
ter of the evolution is also confirmed by the shape of the 
Probability Density Function (PDF), P(T(*)) , which has 
a Poisson-type tail but is also not much different from an 
inverse Gaussian - a typical distribution encountered in 
neural dynamics . 

Finally, given the evidence of the effective stochastic 
behaviour, we develop a standard mean-field approach 
to characterize the behaviour of the network in the large 
iV-limit. In particular, one can conveniently approximate 
the mean field with a deterministic drift plus a zero- 
average stochastic process. This amounts to replacing 
the initial model with the following stochastic equation 
(to be interpreted in the Ito sense, as the noise arises 
from a discrete-time equation), 

G(v.+w) G(vj+w) ,,, ^, 

(22) 

where the time T is determined self-consistently, while 
'rij{t) is a delta-correlated stochastic process 

with Sij and 6(t) being the Kroenecker and Dirac 6- 
functions, respectively. In the thermodynamic limit the 
noise term vanishes and can therefore be neglected in 
the computation of average properties. In particular, by 




FIG. 11: Patterns associated with the two different phases 
of the dynamical system. On the vertical axis (from top to 
down) is reported the time, while the horizontal axis repre- 
sents the index of successive firing neurons. The upper pat- 
tern refers to LP for G = 0.5 , and the lower to the transition 
region with G = 1.3. In both cases a system with = 801 
and dilution r,n = 0.05 is considered. 

integrating the deterministic force field in Eq. 122|l . one 
obtains an implicit equation for T{G), 

„ a; , c — xw 

^ = ; ^ j 

l + x c— 1 — a:;(l+w) 

where x = G/T{G). This formula holds not only in the 
UP, but also in the LP, because the main difference be- 
tween the two regimes does not concern the noise ampli- 
tude, but its correlation properties. In fact, the relative 
diffusion of the generic neurons i and j depends on the 
fact that the stochastic signals ?/; and r]j are not corre- 
lated with one another. If this is the case, the relative 
diffusion induced by the noise over a time equal to the 
average ISI is 1/Vn. Since, on the other hand, the av- 
erage distance between consecutive neuron potentials is 
on the order of 1/A^, we conclude that for any finite G 
and large enough N the diffusion is sufficiently strong 
as to scramble the neurons. The smallness of the noise 
amplitude also explains why at a coarse-grained level the 
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FIG. 12: The normalized autocorrelation Cisi of the T**^ dur- 
ing transient for a single neuron i vs. the integer delay for 
N = 1000, G = 2.0 . The inset shows the corresponding dis- 
tribution of the T^'' and an inverse Gaussian with the same 
mean and variance. 



firing patterns appear almost periodic (see e.g. Fig. |SJ). 
Now, in order to complete a self-consistency argument, 
it would be necessary to connect the changes in the neu- 
ron ordering with the correlation properties of the noise 
terms in Eq. (|22|l . Since different neurons have differ- 
ent connection trees, r/i [t -\- T) will be decorrelated from 
rjiit) and, more important, rii(t + T) — rij{t + T) will 
change with respect to the previous ISI. However, trans- 
forming these hand- waving arguments into a quantitative 
criterion for the identification of the critical value Gc, 
above which self-generated fluctuations can be robustly 
sustained, is not an easy task and we leave it to future 
studies. Here, we limit ourselves to stressing the relevant 
role played by correlations between neighbouring neurons 
(i.e. those with the closest -to-threshold membrane po- 
tentials). Below Gc the neurons rearrangement continues 
until it stops when a suitable ordering is reached, char- 
acterized by the property that neighbouring neurons do 
not diffuse away from each other. The existence of expo- 
nentially long transients in the UP tells us that such an 
arrangement does not exist in the strong coupling regime. 

The two-phase scenario here above described persists 
for a wide range of values. While no qualitative 
changes are expected to appear until becomes so large 
that the network is decomposed into uncoupled blocks, 
we are unable to formulate any conjecture for the limit 
case Tm — > 0, due to the need of studying very large 
networks. 



VI. CONCLUSIONS AND OPEN PROBLEMS 

The main result of this paper is the observation that 
a network of leaky integrate-and-fire neurons can exhibit 
a pseudo-chaotic behaviour in spite of an entirely nega- 



tive Lyapunov spectrum. Previous studies|l7'| have iden- 
tified in strong localized nonlinearities a necessary (and 
far from sufficient) condition for this phenomenon to ex- 
ist. Discontinuities are indeed responsible for a sudden 
amplification of (finite-amplitude) perturbations and for 
the resulting irregular dynamics. In the models where 
this phenomenon was initially observed, the discontinu- 
ities are somehow artificial features of the local maps, but 
it has been recognized that they may also spontaneously 
emerge. This is the case of the hard-point diatomic gas, 
where discontinuities arise in the vicinity of three-body 
collisions [l^. Here, we have seen that they also occur 
in neuronal networks, in connection with changes in the 
spike ordering. However, when and why a discontinuity 
may be so important as to steadily sustain (in the ther- 
modynamic limit) an irregular dynamics is still a com- 
pletely open problem. 

Another open question concerns the nature of the 
phase-transition. Macroscopically, it manifests itself as 
a collective de-synchronization, but there are difficulties 
in quantitatively describing the phenomenon. On one 
hand, it cannot be analysed by looking at the linear sta- 
bility of some solution, since linear stability is ensured 
above and below the transition. Moreover, the absence 
of local instabilities (standard deterministic chaos) makes 
it questionable the concept of a transition point. In fact, 
the study of a partially analogous transition in a coupled 
map lattice revealed that regular and irregular phases 
are separated by a finite fuzzy region, rather than by a 
point-like transition |23| • Perhaps the almost global cou- 
pling of the network studied in this paper may induce a 
standard transition scenario and the analogies with con- 
tact processes suggested by Fig. |21 indicate that this is 
a reasonable expectation. However, for the moment, we 
must limit ourselves to observing that the scenario is ro- 
bust against various modifications, such as the addition 
of further disorder accounting for a non homogeneous 
strength of the coupling and specific modifications of the 
network topology. 

Recently, it has been pointed out that an irregular dy- 
namics should significantly enhance information process- 
ing |23|. At variance with sparsely coupled chaotic net- 
works, where chaos results from the balance of inhibition 
and excitation p^ . we have found irregular firing pat- 
terns, characterized by a Poisson-like distribution of in- 
terspike intervals, in the presence of linear stability and 
absence of excitation. It would be worth investigating 
if and how such a Lyapunov stable irregular dynamics is 
an effective tool for information processing. In the litera- 
ture there is a growing evidence of recurrent motifs in the 
firing pattern (see e.g. i23j), which suggests that phase 
lockings may be a tool to encode information. In our 
model, recurrent motifs appear naturally in the UP and 
yet keep evolving, indicating a sort of information pro- 
cessing which certainly deserves a thorough investigation 
from this point of view. 
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